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In this paper we report on numerical studies of the Cauchy problem for equivariant wave 
maps from 2 + 1 dimensional Minkowski spacetime into the two-sphere. Our results provide 
strong evidence for the conjecture that large energy initial data develop singularities in 
finite time and that singularity formation has the universal form of adiabatic shrinking of 
the degree-one harmonic map from into . 



1 Introduction 

A wave map is a function from the Minkowski spacetime (]R."+^, r;) into a complete Riemannian 
manifold {N,g), U : M"+^ N, which is a critical point of the action 



The associated Euler-Lagrange equations constitute the system of semilinear wave equations 



where F's are the ChristofFel symbols of the target metric g. Wave maps are interesting both 
for mathematicians by providing a simple geometric setting for studying the problems of global 
existence and formation of singularities and for physicists (who call them sigma models) as 
toy models of extended structures in field theory (see for a recent review). In this paper 
we consider the case where the domain manifold is 2 -f 1 dimensional Minkowski spacetime, 
M = R^+^j and the target manifold is the 2-sphere, N — S'^, with the standard metric 
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Abstract 




(1) 



nu^ + Tic{U)daU^d''U^ = 0, 



(2) 



(3) 



We restrict attention to equivariant maps of the form 



u = u{t, r), 9 — (j), 



(4) 



where (r, are the polar coordinates on M^. The wave map system (|]) reduces then to the 
semihnear scalar wave equation 



1 

Utt = Urr H U,. — 

r 



sin(2M) 



(5) 
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The main open question for equation (|5|) is the issue of global regularity, namely: do all solutions 
starting with smooth initial data 



remain smooth for all times, or do they lose regularity for some data? Our paper reports on 
numerical investigations of this problem. 

Note that the conserved energy associated with solutions of (0) 



is invariant under dilations: if u\(t,r) — uit/ X,r / X), then i?[uA] = E[u\. In this sense 2+1 is a 
critical dimension for wave maps. Note also that the requirement that energy be finite imposes 
a boundary condition at spatial infinity oo) — kir {k — . . .) which compactifies into 

and thus breaks the Cauchy problem into infinitely many disconnected topological sectors 
labeled by the degree k of the map ^ S^. 

Let us recall what is known rigorously about the problem. Besides the routine local existence 
proof, the only global result is that there is a unique smooth solution for all times provided 
that the initial energy is sufficiently small 0|. In order to obtain global existence without the 
assumption of small energy it would be sufficient to show that the energy cannot concentrate at 
the hypothetical singularity. This kind of result was derived by Shatah and Tahvildar-Zadeh 
for convex targets and then extended by Grillakis |^ for nonconvex targets with arbitrarily 
narrow neck. However, the Morawetz identity, which is the basic tool in proving nonconcentration 
of energy, loses the desired positivity properties for the geometry of two-sphere. This raises the 
question: is the restriction on admissible targets assumed in |3| only of technical nature or is it 
essential? In particular, is concentration of energy possible in the case of two-sphere as target? 

We present numerical evidence that for large energy solutions the energy does in fact concen- 
trate and consequently the solutions blow up in finite time. We show that the process of energy 
concentration proceeds via adiabatic evolution along the one-parameter family of dilations of the 
degree-one static solution. In this sense the shape of blowup is universal. The rate of blowup, 
determined by the asymptotic speed of adiabatic evolution, is slower than that predicted by the 
geodesic approximation and goes to zero as the singularity is approached. 

The paper is organized as follows. In section 2 we derive static solutions and study their 
stability. As mentioned above these solutions play an essential role in the process of singularity 
formation. For completeness, in section 3 we discuss singular self-similar solutions and argue 
that they play no role in the Cauchy problem. The main body of the paper is contained in 
section 4 where we present the results of numerical investigations. On the basis of these results 
we formulate three conjectures about the nature of singularity formation in the model. Finally, in 
section 5 we comment on some earlier work on this problem and point out some open problems. 



u{0,r) = ua{r), Ut{0,r) ^ Ui{r) 
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2 Static solutions 



Static solutions of equation (j^) can be interpreted as spherically symmetric harmonic maps from 
the Euclidean space into S^. They satisfy the ordinary differential equation 

,, 1 , sin(2M) , . 

+ - - 0, (8) 

where ' — d/dr. The obvious constant solutions of (||) are u — and w = tt; geometrically these 
are maps into the north and the south pole of 5*^, respectively. The energy of these maps attains 
the global minimum E — 0. Another constant solution is the equator map u — 7r/2 but this 
solution is singular and has infinite energy. The fact that equation is scale invariant does 
not exclude nontrivial regular solutions with finite energy (Derrick's argument is not applicable) 
and, in fact, such solutions are well-known both in the mathematical literature as harmonic maps 
from into S'^ and in the physics literature as instantons in the two-dimensional euclidean 
sigma model. They can be derived in many ways. One possibility is to use a Bogomol'nyi-type 
argument which goes as follows. Let x = Inr and U{x) = u{r). Then, assuming staticity 

oc 2 oo oo 

E[u] =71 J {u'^ + rdr = n J {U'^ + sin^C/) dx = tt J (U' - sinU f dx - 27rcos U 

— oo — oo 

(9) 

Thus, in the topological sector fc = 1 the energy attains the minimum, E = An, on the solution 
of the first order equation U' — sinC/ = 0, which is U{x) = 2arctan(e^) up to translations in x. 
Therefore 

us{r) = 2arctan(r) (10) 

is the static degree-one solution (the problem has reflection symmetry so, of course, —us{r) is 
also the solution). By dilation symmetry, the solution us{r) generates the orbit of static solutions 
Ug{r) = us{r/X). We remark in passing that the solution ( |l0[ ) can be alternatively obtained 
in the elegant geometric way by taking the identity map between two-spheres and making the 
stereographic projection. 

We consider now the linear stability of the static solution us{r). Inserting u{t,r) = us(r) + 
e'*'*f (r) into (|^) and linearizing, we get the eigenvalue problem (the radial Schrodinger equation) 

Lv=(-^--^ + V{r))v = ev, (11) 
dr'' r dr 



where 



cos(2w5) 1 — 6r^ 

(1 -t- r2)2r2 



This potential has no bound states as can be shown by the following standard argument. Consider 
the perturbation induced by dilation 



d 



dX 



2r 



^oW = "TT^sW , , = ru'sir) = (13) 



A=i ' 1 + r 
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This is the solution to = (so called zero mode). The fact that the zero mode vo{r) has 
no nodes implies by the standard result from Sturm-Liouville theory that there are no negative 
eigenvalues, and eo ipso no unstable modes around us{r). Note that the zero mode is not a 
genuine eigenfunction because it is not square integrable. Therefore the operator L has the 
purely continuous spectrum fc^ > 0. 



3 Nonexistence of self-similar solutions 

Since equation (jsl) is scale invariant, it is natural to look for self-similar solutions of the form 



«(i,0 = /(^^j, (14) 

where T is a positive constant. Substituting this ansatz into ^ one obtains the ordinary 
differential equation 

, A P \df sin(2/) ^^^^ 



dp^ \p J dp 2p2(l -p2) 



Let us consider equation (|1^) inside the past light cone of the point (t = T,r = 0), that is for 
p G [0, 1]. It is well-known that there are no solutions which are analytic at the both ends of this 
interval [Q. However, it is not-well-known that there arc solutions which are less regular. They 
can be easily derived by setting p = 1/ cosh(?/) , so that (|l^) simplifies to 



dy 



2 sin(2/) = 0. (16) 



This equation is solved (up to translations in y) by f{y) = 2arctan(e^), so in terms of p we get 
a one-parameter family of self-similar solutions 

Up) - 2arctan ( ) • (17) 

V1 + v/1-pV 

These solutions are analytic at p = but they are not differentiable at p = 1 (and consequently 
have infinite energy). Since such solutions cannot develop from smooth initial data inside the 
whole past light cone of the singularity, they are not expected to play any role in the Cauchy 
problem. The numerical results described below support this expectation. In this respect the 
wave maps in 2-1-1 dimensions are completely different from the wave maps in 3-1-1 dimensions 
where a stable analytic self-similar solution determines the process of singularity formation M . 



4 Numerical results 

We have solved numerically the Cauchy problem (|^), ^ for various one-parameter families of 
initial data which interpolate between small and large energy. The details of our numerical 
methods are described in the appendix. The results described below are universal in the sense 
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that they do not depend on the choice of the family of initial data (nor on the topological sector) . 
For concreteness we present them for the degree-zero initial data of the form 



^(0,r) = a(-^)'cxp -(^^) 



ut(0,r)=0, (18) 



where the amplitude A, the radius i?, and the width S are free parameters. Below we fix i? = 2, 
S = 0.4 and vary A. We emphasize that the amplitude is by no means distinguished: any 
parameter which controls the energy of the initial data could be varied. Note that regularity at 
the center requires that u{t, 0) = for all t < T, where T is a time when the first singularity 
(if there is any) develops at r = 0. Since the initial momentum is zero, the initial profile splits 
into two waves, ingoing and outgoing, traveling with approximately unit speed. The evolution 
of the outgoing wave has nothing to do with singularity formation so we shall ignore it in what 
follows. The behavior of the ingoing wave depends on the amplitude A. For small amplitudes, 
the ingoing wave shrinks, reaches a minimal radius, and then expands to infinity leaving behind 
the zero energy region (see figure 1). When A increases, the minimal radius at which the wave 
bounces back decreases and seems to go to zero for some critical value of the amplitude A*. 
Finally, for the supercritical amplitudes A > A* , the wave does not bounce back and keeps 
shrinking to zero size in finite time. More precisely, we observe that evolution of the wave near 
the center (so called inner solution) is well approximated by the degree-one static solution us 
with a time-dependent scale factor A 

u{t, r) « -2 arctan (^1^^ • (19) 

We shall refer to this formula as the adiabatic approximation. Using the adiabatic approximation 
the evolution of the ingoing wave can be described as follows. For subcritical amplitudes A < A* 
the scale factor X{t) decreases, attains a minimum Amm, and then increases. When A A* , 
then Xmin 0. For supercritical amplitudes A> A* , the scale factor decreases monotonically to 
zero in finite time. As follows from (|l9|), Mr(i, 0) ~ A~^(t), hence for the supercritical solutions 
the gradient blows up at the center in finite time. Various aspects of this behavior and some 
numerical details are shown in figures 2, 3, 4, and 5. 

We have not been able to develop a rigorous mathematical understanding of the adiabatic 



approximation (19) but we can make some heuristic arguments which help understand the ob- 
served behavior. Let us define a similarity variable 77 = r/\[t). Substituting u = u{t,r]) into 
equation (||) gives 

-ii + (1 - \W)u" + [l + (AA - 2X')r^'] - - = 0, (20) 

where ' = d/dt, ' = d/drj. In order to "explain" the observed behavior we make two assumptions. 
The first assumption, which is the essence of adiabaticity, says that the dynamics of the solution 
near the center is slaved in the varying scale A(t). This implies that we can neglect the term ii 
in (^o|) which refers to the explicit time dependence. The second assumption concerns the rate 
of blowup and says that 

^0 as t /T. (21) 

T-t ^ ^ ' 
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Eq.(pTl) implies that the terms involving the time derivatives of A in Eq.(pO|) tend to zero as 
t T . After dropping these terms and the u term, Eq.(20) becomes formally the same as 
Eq.(||) and therefore is solved by us(j]). Although this ad hoc explanation is certainly not 
satisfactory, it is fully consistent with the numerics. In particular, it explains why the "amount" 
of the inner solution which is approximated by ( ]l9| ) increases as the wave shrinks (see figure 5). 

Of course, it would be very interesting to find an exact rate of blowup. We obtain a reasonable 
fit to the power law behavior 

A(t) - (T - t)" as t /■ r, (22) 
with the anomalous exponent a « 1.1 ± 0.05 (see figure 6). However, in view of the limited 



resolution of our numerics near the blowup and the lack of theoretical arguments behind (22) 
we caution the reader not to take ( ^^ as a serious prediction; in particular we cannot rule out 
logarithmic corrections to the power law behavior. 

Now, on the basis of the numerical studies just described we would like to put forward three 
conjectures which summarize the main points of our findings. 

Conjecture 1 (On blowup for large data). For initial data with sufficiently large en- 
ergy, the solutions of equation (^ blow up in finite time in the sense that the derivative Ur{t, 0) 
diverges as t T for some T > 0. 

Conjecture 2 (On blowup profile). Suppose that the solution u{t, r) of the initial value prob- 
lem blows up at some time T > 0. Then, there exists a positive function X(t) \ for 
t y T such that 

lim u{t, Xit)r) = ±us(r) for r > 0. (23) 

Conjecture 3 (On energy concentration). Suppose that the solution u{t,r) of the initial 
value problem ^,^) blows up at some time T > 0. Define the kinetic and the potential energies 
at time t < T inside the past light cone of the singularity by 

T-t T-t 

/r sin u 
ulrdr and Ep{t) = n (ul -\ —) rdr. (24) 



Then: 

(i) the kinetic energy tends to zero at the singularity 

Yon E Kit) ^0; (25) 

t/^T 

(ii) the potential energy equal to the energy of the static solution ug concentrates at the singu- 

larity 

lim Ep(t) = E\us\ = 4tt. (26) 
tyT 



6 



We have already discussed the evidence for Conjectures 1 and 2. Conjecture 3 is basicaUy the 
consequence of (|l|) and @. To see this let us substitute us{r/X(t)) into ( |2^ ) to get 



EK{t)^TTX^ j u's^r^dr and Ep{t) = tt j {u's^ + ^^^^) rdr. (27) 



Assuming (pl|), the upper limits in these integrals tend to infinity as t /" T, so ( p5|) and (26) 
follow (notice that the integral in Ek diverges logarithmically). Conjecture 3 means that as the 
blowup is approached the excess energy above the energy of the static solution flows outward 
from the inner region. This is clearly seen in our simulations. 
We remark that an averaged weak version of (Eq) 



lim i J EK{t)dt = (28) 

T-c 

was proved by Shatah and Tahviladar-Zadeh Q . 



5 Final remarks 

We would like to comment on two papers by Piette and Zakrzewski |^ and by Linhart |^] which 
were devoted to the adiabatic evolution in the degree-one topological sector. These authors 
considered initial data which have the shape of the static solution and a nonzero momentum 
directed inwards. They observed adiabatic shrinking (^^ with the scale factor A changing 
approximately linearly in time. This was "explained" by the geodesic approximation (an old 
idea due to Manton [||) as follows. Substituting the ansatz (|l^) into the action one obtains the 
effective action for the scale factor \{t). The potential energy part does not depend on A, so 
only the kinetic energy part contributes to the effective action, 

oo oc ^ 

S,ff[\] = j dt j uirdr^ j dt }? J ~ c J X^dt, (29) 



where the "constant" c is logarithmically divergent^. As long as the geodesic approximation is 
used to model blowup, this divergence is irrelevant because it can be removed by truncating the 
action at some large radius. By finite speed of propagation, such a truncation cannot affect the 
blowup. Thus, the effective action (|2^) yields the scale factor going to zero linearly A ~ T — i 
as i y T. To verify the accuracy of this approximation, Piette and Zakrzewski have solved the 
Cauchy problem numerically and got X {T — t)^+'^, where e is a positive number of the order 
of 0.1. The smallness of e was interpreted in favor of the geodesic approximation. Although 
we have confirmed these results numerically (see the discussion above and figure 6), we disagree 

^In the language of geodesic approximation this divergence means that the volume of the moduli space is 
infinite which nota bene is equivalent to the fact that the zero mode is not square integrable. In the literature 
one can find statements that zero modes which are not square integrable are "frozen" by infinite inertia. The 
papers ^ and demonstrate that these statements, based on the naive picture of ( ps| ) as the action for the free 
particle with infinite mass, are wrong. 
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with the authors of ||] and Jzt regarding the accuracy of the geodesic approximation. As we 
wrote above, the assumption (21) is crucial for the vahdity of the adiabatic approximation. In 
contrast, the hnear decay of A, predicted by the geodesic approximation, is inconsistent with the 
observed adiabatic evolution, and, as follows from (p7|), gives the wrong prediction on energy 
concentration at the blowup which contradicts Conjecture 3 and theorem (28). The reason for 
which the geodesic approximation fails to capture these crucial features of blowup is easy to 
understand: this approximation completely neglects radiation which is essential in expelling the 
excess energy from the inner region. 

Besides the obvious problem of proving Conjectures 1, 2, and 3, the research presented here 
raises a number of questions. Most interesting among them, in our opinion, are: 

1. What mechanism selects the time evolution of the scale factor A? In answering this question 
the methods of center manifold theory might be useful. In particular, using weakly nonlinear 
stability analysis it should be possible to derive the amplitude equations for the nonlinear evo- 
lution of the zero mode. The problem shares many features with the problem of blowup for the 
nonlinear Schrodinger equation in two spatial dimension. It is feasible that the techniques of 
asymptotic matching used there [|o) could be also applied to our problem. 

2. What is the evolution at the threshold for singularity formation? What does the fine-tuning 
accomplish dynamically? The model does not fit into the framework of "standard" threshold 
behavior where a codimension-one stable manifold of a certain critical solution separates blowup 
from dispersion Unless more accurate fine-tuning would reveal a new universal behavior 
very close to the threshold (which is unlikely), it seems that solutions evolving from exactly 
critical initial data also blow up in the adiabatic manner but at a much slower rate (see Fig. 6). 

3. To what extent the results are specific to the equivariance ansatz; is the blowup stable under 
general perturbations? 

We hope to be able to say more about these problems in future publications. 
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Appendix: Numerical methods 

In order to solve equation (^ numerically we rewrite it as the first order system in time: 

ut = V, (A. la) 

1 sin(2M) / A 11 ^ 

Vt = Urr+-Ur ^ „ ■ (A. lb) 

r Ir'^ 

We solve this system by finite differencing. To ensure regularity at the origin we require that 
u(r, t) — 0{r) for r — > 0, from which the inner boundary conditions follow: u(0, t) ~ 0, w(0, t) — 0. 
As the outer boundary condition we impose an approximate outgoing wave condition. A naive 
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centered difference scheme applied to the right hand side of (A. lb) would trigger an instability 
near r = 0. To avoid this, we use a scheme which is natural for the operator 



Cu = Urr + -Ur = -dr{rdrU), (^-.2) 



and takes the form 

Cu « - 
r 



1 A h u{r + h) — u{r) h u{r) — u{r — h) 



hV 2' h '2' h 



(A.3) 



where /), denotes the spatial mesh size. 

For time evolution we use a standard leapfrog scheme. In that way we obtain a scheme which 
is second order accurate in space and time. 

Preliminary results obtained on a imiform, fixed grid show that the most interesting inner 
solution is well approximated by the static solution us with a time-dependent scale factor A 
(see equation (19)). Therefore, to follow this solution we have to change the scale of spatial 
resolution in time and keep it roughly proportional to A, at least as long as A decreases. To this 
end we apply an adaptive algorithm in which both the mesh size and the time step are refined as 
the solution shrinks. We start with a uniform grid covering an interval (0, R) with some initial 
resolution characterized by a mesh size h = Aro. We continue the evolution of the system on 
this grid with a time step A^o as long as Ur{0,t)h < C where C is some small fixed constant 
(spatial tolerance factor). When this inequality is violated we refine the original grid on the 
interval (0, R/2) by covering it with the resolution Ari = Aro/2. The values of functions at the 
new points, not defined on the parent grid, are obtained by interpolation. From that time on we 
continue the evolution on the finer grid with the time step Ati = Ato/2. Iterating this process 
several times wc get the resolution adaptively adjusted to the solution. 

In order to make sure that the numerical results are reliable, we have reproduced them using 
a different implicit finite-differencing scheme in which we used ln(r) as the spatial variable. 
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Figure 1: Snapshots of the evolution of initial data (|^) with small amplitude A — 0.5. The 
ingoing wave bounces back and disperses. The minimal radius is attained at t w 2.4. 
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Figure 2: The scale factor X{t) for a sequence of solutions with nearly critical initial amplitudes. 
Numerically, X{t) was calculated from the formula Ur{t,0) = —2/\{t). The critical amplitude 
was estimated to be A* w 1.0678281. 




_2 g I , , u , , U , , , , , ,__._L 

ie-007 1e-006 1e-005 0.0001 0.001 0.01 0.1 1 

r 




_2 g I , , u , , u , . u , , u , , u , , ,_ 

ie-007 1e-006 1e-005 0.0001 0.001 0.01 0.1 

r 

Figure 3: Snapshots of the evolution for the marginahy subcritical amphtude A = 1.06782. 
The upper plot shows the shrinking phase while the lower plot shows the expanding phase. 
The horizontal arrows indicate the direction of motion. For sufficiently small r all profiles 
have the shape of the suitably rescaled static solution. The scale factor attains the minimum 
Amm « 0.9308 X 10^5 for t « 2.6086. Notice that u{t, r) > -tt for aU times. 
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Figure 4: The momenta for the same data as m Figure 3, rescaled by the factors X{t). In 
agreement with the adiabatic approximation ( p^ the profiles have the shape of the zero mode 
vo{r) with the amphtudes given by the logarithmic derivative of the scale factor. 
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Figure 5: Snapshots of the supercritical evolution for the amplitude A — 1.072. The upper plot 
shows the ingoing wave shrinking indefinitely as t approaches the blowup time T « 2.5593558. In 
the lower plot the profiles from the upper plot, rescaled by the factors X{t), are shown to collapse 
to the static profile —us{r). Notice that the fifth profile, corresponding to T — t = 0.0002888, 
overshoots — tt; this seems to be the necessary and sufficient condition for blowup, which is 
reminiscent of a similar phenomenon for the heat flow for harmonic maps 1^ . 
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Figure 6: The scale factor from two supercritical evolutions. In both cases the asymptotic 
behavior of X{t) is well approximated by the power law A (T — f)" with the exponent aw 1.1. 
However, in the marginally supercritical case the transient regime with a w 2.3 in clearly seen 
before the true asymptotic state is reached. As A ^ A*, the crossover between the transient 
and the asymptotic regimes occurs closer and closer to the blowup, which suggests that solutions 
with exactly critical initial data blow up at the much slower "transient-regime" rate. 



